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Abstract: The genus Trachypithecus is the most diverse langur taxon, distributed in southwestern China, south and southeastern Asia. 


In this study, we include 16 recognized Trachypithecus species to reconstruct the phylogeny with particular concern to the taxonomy 


of the three subspecies of T. phayrei using multiple genes. Our results support a sister-relationship between T. p. phayrei and T. p. 


shanicus. However, the mitochondrial CYT B gene supported T. p. crepuscula as a distinct species, but the nuclear PRM/ gene 


suggested a closer relationship between T. p. crepuscula and T. p. phayrei. The incongruence between nuclear and mitochondrial 


genes suggests that hybridization may have occurred, a fact that would benefit from re-examination using multiple unlinked nuclear 


genes. 
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The genus Trachypithecus is the most diverse langur 
taxon, having a broad distribution including India, Sri 
Lanka, Bangladesh, Southwestern China, and Southeast 
Asia (Groves, 2005; Wang et al, 1999). It is 
phylogenetically embedded within the Family 
Cercopithecidae and closely related to Semnopithecus 
(Perelman et al, 2011; Wang et al, 2012). Groves (2005) 
assigned full species status to 17 taxa, which he clustered 
into 5 species groups. While 16 of these species have 
been assessed in other phylogenetic contexts (Bleisch et 
al, 2008; Geissmann et al, 2004; Karanth et al, 2008; 
Liedigk et al, 2009; Nadler et al, 2003; Wang et al, 2012; 
Zhang & Shi, 1993), which has greatly improved our 
understanding of Trachypithecus evolution, there has 
been no dedicated molecular analysis of the complex 
species structure of Trachypithecus. Of particular 
concern is the disputed taxonomy of the endangered 
Phayre’s leaf monkey (T. phayrei). Three putative 
subspecies inhabit Bangladesh, northeastern India, 
Myanmar, Southwestern China, Thailand, Laos, and 
northern Vietnam (Bleisch et al, 2008). Recent genetic 
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analyses have demonstrated that T. p. crepuscula and T. p. 
phayrei do not form a monophyletic clade. T. p. phayrei 
from India is the sister taxon of T. barbei and T. obscures, 
but T. p. crepuscula from Vietnam represents a distinct 
lineage, being a close relative of the T. francoisi species 
group (Karanth et al, 2008; Nadler et al, 2003). 
Accordingly, Nadler et al (2003) and Liedigk et al (2009) 
suggested that T. p. crepuscula from Vietnam should be 
given full species status. However, a lack of genetic 
information from T. p. crepuscula and T. p. shanicus 
(southwestern China) has prevented agreement on the 
purported taxonomy (Bleisch et al, 2008). 

In this study, we sampled T. p. crepuscula and T. p. 
shanicus from southwestern China (Yunnan) and 
northern Myanmar and sequenced both mitochondrial 
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and nuclear genes to resolve the phylogenetic 
relationships of T. phayrei. The sequences of T. p. 
shanicus are represented for the first time. By including 
langur sequences from GenBank into our analysis, we 
test the validity of the species groups proposed by 
Groves (2005). Furthermore, we estimated divergence 
times simultaneously with phylogenetic reconstruction to 
test the diversification pattern of Trachypithecus through 
time (Drummond et al, 2012). This effort will provide 
further evidence to clarify the nebulous classification of 
T. phayrei and the putative Trachypithecus species-group 
divisions. 


MATERIALS AND METHODS 


Sampling and lab work 
Four samples of T. p. crepuscula and five of T. p. 
shanicus were collected from northern Myanmar, Yunnan, 


and Guangxi, China. One sample of T. germaini was 
obtained from the Kunming Cell Bank of the Chinese 
Academy of Sciences (Table 1). These samples are 
tissues from deceased individuals, museum skins, or 
noninvasive feces collections. All animal samples were 
obtained following the regulations of China for the 
implementation of the protection of terrestrial wild 
animals (State Council Decree [1992] No. 13) and 
approved by the Ethics Committee of Kunming Institute 
of Zoology, Chinese Academy of Sciences, China. 

Total DNA was extracted from tissues and museum 
skins using the phenol/proteinase K/sodium dodecyl 
sulphate method (Sambrook & Russell, 2001). Four fecal 
samples were collected from the Gaoligong Mountain 
National Nature Reserve in 2008. Samples were stored 
using the “two-step” storage procedure (Nsubuga et al, 
2004). We extracted total DNA from fecal samples using 
the 2CTAB/PCI method (Vallet et al, 2008). 


Table 1 Sample and sequences used in this study 








Species Collection code Code Sample Locality CYTB PRMI1 LZM 

p ae EUA la Myanmar01 NM Northern Myanmar KC285863 KC285882 KC285875 

YN_WLS0301001 WL Mt. Wuliang, Yunnan,China KC285866 KC285883 KC285876 

YN_XSBNI1 BNI Xishuangbanna, Yunnan,China KC285864 KC285884 KC285873 

GX_KCB89 BN2 Xishuangbanna, Yunnan,China KC285865 KC285885 KC285874 

T. p. shanicus YN_GLG0912002 G1 Mt.Gaoligong, Yunnan,China KC285867 KC285886 KC285877 

YN_GLG08002 G2 Mt.Gaoligong, Yunnan,China KC285868 KC285887 

YN_GLG08003 G3 Mt.Gaoligong, Yunnan,China KC285869 = KC285878 

YN_GLG08004 G4 Mt.Gaoligong, Yunnan,China KC285870 = KC285879 

YN_GLG08005 G5 Mt.Gaoligong, Yunnan,China KC285871 = KC285880 

T. germaini KCB92007 V2 Northern Vietnam KC285872 KC285888 KC285881 
We sequenced the mitochondrial CYT B and nuclear Inc., USA) and aligned with MUSCLE (Edgar, 2004). 


protamine P1 (PRMI) and lysozyme (LZM) genes in this 
study. CYT B has been widely used in mammalian 
phylogenetic and phylogeographic studies (Bradley & 
Baker, 2001). This gene was amplified using L14724_hk6 
(CCGTGATATGAAAAACCATCGTTG) and H15915 
(Irwin et al, 1991). L14724_hk6 was modified from the 
universal primer L14724 (Irwin et al, 1991) to avoid 
amplification of nuclear pseudogenes of CYT B (Karanth, 
2008). PRM1 and LZM have been used previously in the 
phylogenetic study of other primates, including langurs 
(Karanth et al, 2008). These two genes were amplified 
following (Karanth et al, 2008). All PCR products were 
purified using the UNIQ-10 spin column DNA gel 
extraction kit (Sangon, Shanghai, China). Purified products 
were directly sequenced with PCR primers using the 
BigDye Terminator Cycle kit v3.1 on an ABI 3730xl 
sequencer by Tiangen Biotech Co, LTD. (Beijing, China). 


Phylogenetic and molecular dating analyses 
Nucleotide sequences were edited using SeqMan 
and EditSeq in the DNASTAR package v7.1 (DNASTAR, 
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CYT B and coding regions of the two nuclear genes were 
translated into amino acids following the identification of 
any premature stop codon. All these sequences were 
submitted to GenBank (Accession numbers: KC285863 
- KCKC285888). 1141 bp of CYT B, 387 bp of PRM/, 
and 135 bp of LZM were used in further analyses. 
Additional sequences of Catarrhini were downloaded 
from GenBank and added to the alignments 
(Supplementary Table 1, available online). In total, 70 
CYT B, 35 PRM1/ and 32 LZM sequences representing 36 
species including all three subspecies of Trachypithecus 
phayrei were used in this study. 

We performed partitioned Bayesian analyses 
(Brandley et al, 2005) to reconstruct the phylogenetic 
relationships using BEAST v1.7.2 (Drummond et al, 
2012). The analyses were performed on both the CYT B 
data set as well as the three-gene combined data set. 
Missing data were coded as “?”. The model of DNA 
evolution was determined by a Bayesian Information 
Criterion (BIC) (Schwarz, 1978) in jModelTest v0.1.1 
(Guindon & Gascuel, 2003; Posada, 2008). Because 
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there are missing data in the CYT B alignment, this gene 
was divided into two partitions (partition 1: 1-423 bp, 
826-1 141 bp; partition 2: 424-825 bp). Because 
partition 2 contains no missing data, it was used for the 
estimation of the time of molecular divergence (see 
below). The evolutionary model for each codon position 
of partition 1 and partition 2 of CYT B, LZM, and the 
coding and non-coding regions of PRMI/ were 


CYTB 










determined separately. We did not calculate the 
evolutionary model for each codon position of LZM or 
the coding regions of PRMI to avoid error caused by 
overpartitioning (Brown & Lemmon, 2007). In 
jModelTest, three substitution schemes were selected and 
a proportion of invariant sites were not included in the 
model selection, following He et al (2012). The partition 
information and models selected are shown in Figure 1. 


1st codon: HKY+G 
2nd codon: HKY+G 
3rd codon: GTR+G 


1-423 424-825* 826-1141 


Ist codon: GTR+G 





2nd codon: HKY+G 
3rd codon: HKY+G 


PRMI 








coding: HKY 


non-coding: SYM+G 


LZM 


K80+G 


Figure 1 Models of molecular evolution used for each gene/partition for phylogenetic reconstruction 


*: partition was used for divergence time estimation. 


We performed a MCMC search of twenty million 
generations, sampling every 2 000 generations. We 
constrained the monophyly of Semnopithecus + 
Trachypithecus according to previous analyses 
(Osterholz et al, 2008; Perelman et al, 2011; Wang et al, 
2012). All analyses were repeated twice. Tracer v1.5 was 
used to make sure all analyses reached the same posterior 
distributions and estimated the convergences by 
calculating effective sample sizes (ESSs) (Rambaut & 
Drummond, 2009). Posterior probabilities (PP) > 0.95 
were considered as strongly supported (Huelsenbeck & 
Rannala, 2004). 

Molecular divergence times were estimated 
simultaneously with the multi-gene phylogenetic 
reconstruction. Because missing data will mislead 
estimates of branch lengths and thus affect the result of 
divergence time estimation (Lemmon et al, 2009), our 
molecular dating analyses were limited to partition 2 of 
the CYT B alignment (402 bp, Figure 1). We used a 
“relaxed” molecular approach (Drummond et al, 2006) 
and used fossil records and secondary calibration for 
dating. Six calibration points were applied as lognormal 
or exponential distributions (Ho, 2007). Saadanius 
hijazensis dated to 29-28 Ma for the hominoid- 
cercopithecoid divergence between 29-28 and 24 Ma 
(Zalmout et al, 2010), so that we set the prior as 
lognormal and offset=24, SD=0.98. The oldest known 
hominoid (Morotopithecus bishopi) is dated to 20.6 Ma 
(Gebo et al, 1997). We set the prior as lognormal. The 
earliest sample age was set to 20.6 Ma (offset=20.6) and 
the older 95% CI to the beginning of the Miocene (24.1 
Ma, SD=0.76). The Homo-Pan divergence occurred 
between 7.2—5.6 Ma (Aiello & Collard, 2001; Senut et al, 
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2001), so we set the prior as lognormal and offset=5.6, 
SD =0.29. The oldest Colobinae was dated to 9.8 Ma 
(Benefit & Pickford, 1986; Nakatsukasa et al, 2010), but 
the molecular divergence estimation indicated a much 
more ancient divergence (Chatterjee et al, 2009), so we 
set the prior as exponential. The earliest possible age was 
set to 9.8 Ma, and the older 95% CI to 20 Ma (mean=3.4). 
The oldest African Colobinae was dated to 6.1 Ma 
(Gilbert et al, 2010). We set the prior as lognormal, the 
earliest possible age to 6.1 Ma, and the 95% CI to the 
beginning of the Late Miocene (11.6 Ma, SD=1.04). The 
most recent common ancestor (MRCA) of 
Semnopithecus and Trachypithecus was estimated to 
have existed at about 4.05 Ma (95% CI=2.93-5.36) 
(Perelman et al, 2011), so we set the prior as lognormal 
(mean=4.03, SD=0.18). Fossils of Trachypithecus have 
been found in Southeast Asia and southern China, most 
of which date to the Pleistocene (Jablonski, 2002). The 
oldest known fossil is T. auratus sangiranensis with an 
age of 1.9+0.5 Ma (Jablonski & Tyler, 1999); however, 
the site stratigraphy remains in doubt (Larick et al, 2000). 
Other fossils are not suitable for calibration, because of 
taxonomic uncertainty. 

The minimum-spanning tree of PRM/ were derived 
from Network v4.5 using the median-joining approach 
(Bandelt et al, 1999). 





RESULTS 


Phylogenetic relationships 

The combined phylogenetic tree is shown in Figure 
2, and the mitochondrial phylogenetic topology is 
identical to that based on multiple genes. Trachypithecus 
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Figure 2 Phylogenetic relationships of Cartarrhini based on combined mitochondrial and nuclear genes and chronogram of 


Cercopithecidae based on partial CYT B gene 


Node numbers above branches indicate mean divergence time; numbers below branches indicate Bayesian posterior probabilities; node bars indicate the 95% CI 


for the clade age. 
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taxa fell into two clades. T. pileatus, T. geei from India, T. 
vetulus, and T. johnii cluster with Semnopithecus (PP=1.0), 
which is consistent with previous studies (Chatterjee et al, 
2009; Karanth et al, 2008; Osterholz et al, 2008; Wang et 
al, 2012), suggesting that T. vetulus and T. johnii should 
be merged into the genus Semnopithecus. The rest of the 
Trachypithecus taxa form a monophyletic clade (PP=1.0) 
with five strongly supported lineages recognized 
(PP=1.0). Our results strongly support that T. p. phayrei 
and T. p. shanicus are sister-taxa (PP=1.0) and are close 
relatives to T. barbei and T. obscures (PP=1.0). T. p. 
crepuscula from southwestern China and Vietnam form a 
distinct lineage (PP=1.0) as the sister taxon to the species 
group T. Obscurus +T. cristatus (PP=1.0). T. geei and T. 
pileatus from Bhutan are strongly supported as the sister 
taxa of T. shortridgei (PP=1.0). 

LZM shows low variation in Trachypithecus and 
Semnopithecus. Only two mutations that lead to amino 
acid substitutions were observed (Table 2). As shown 
previously (Karanth et al, 2008), our results indicate that 
the T. pileatus and T. geei from India, which are 
members of the genus Semnopithecus in the phylogenetic 
tree (Figure 2), have the same LZM amino acid 
sequences as the species in the Trachypithecus clade, 
indicating possible hybridization. 

Nucleotide substitutions and insertions/deletions 
(indels) were found in the PRM1 alignment in both 
coding and noncoding regions. In the Network tree, the T: 


© Semnopithecus clade 
O Trachypithecus clade 


=> 1 amino acid deletion 
— 1 amino acid substitution 
— 1 mutation 







geei from India and T. francoisi share the central 
haplotype, which differs from other haplotypes by two to 
nine mutations (Figure 3). These mutations also lead to 
amino acid substitutions/deletions. T. p. shanicus shares 
the same haplotype with T. germaini and T. obscurus and 
differs from the haplotype shared by T. p. phayrei and T. 
p. crepuscula by one amino acid deletion. 


Table 2 The first ten amino acid sequences of the lysozyme 
protein of presbytini species 





T. cristatus 


Nasalis larvatus M K A L I I L G L V 
Pygathrix nemaeus M A L I I LGL V 
Semnopithecus entellus M K A L T I L G L V 
Trachypithecus johnii M K A L T I L G L V 
T. vetulus M K A L T I L G L V 
T. pileatus (India) M R A L I I L G L V 
T. geei (India) M R A L I I L G L V 
T. francoisi M R A L I I L G L V 
T. phayrei phayrei M R A L I I L G L V 
T. phayrei shanicus M R A L I I L G L V 
T. phayrei crepuscula M R A L I I L G L V 
T. obscurus M R A L I I L G L V 
T. germaini M R A L I I L G L V 
T. cristatus M R A L I I L G L V 

T. pileatus (India) i 
2 |" 
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Figure 3 The median-joining network of the genus Trachypithecus based on partial PRM/ 


Numbers by the branches indicate the number of mutations. 


Divergence time estimation 

The estimates of divergence time including 95% 
credible intervals are shown in Figure 2. The 
Trachypithecus clade diverged from outgroups at about 
5.25 Ma (7.01-3.77 Ma). Diversifications within this 
clade started at about 2.35 Ma (3.43-1.46 Ma) and 
occurred throughout the Pleistocene. T. p. phayrei and T. 
p. shanicus diverged at about 0.37 Ma (0.17-0.63 Ma). 
The MRCA of T. p. crepuscula occurred at about 0.41 
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Ma (0.16-0.77 Ma). 
DISCUSSION 


As has been discussed previously, the conflicting 
phylogenetic information observed in Trachypithecus is 
likely related to hybridization, incomplete lineage sorting, 
and ancestral polymorphism (Karanth et al, 2008; 
Liedigk et al, 2009; Wang et al, 2012). Pleistocene 
glaciations may have caused large scale deforestation 
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throughout southeast Asia, isolating formerly widespread 
species in refugia, and the subsequent recolonization of 
these territories during interglacial periods may have 
allowed for the hybridization of these isolates (Brandon- 
Jones, 1996; Meijaard & Groves, 2006). Consistent with 
the karyotype stability within the subfamily Colobinae 
(O'Brien et al, 2006), hybridization between different 
species of Trachypithecus and between Trachypithecus 
and closely related genera has been observed in captivity 
and could explain the paraphyly of T. pileatus and T. geei 
(Que et al, 2007; Schempp et al, 2008; Wangchuk, 2005; 
Wangchuk et al, 2008). Coupled with hybridization, the 
observed discrepancies between autosomal and 
mitochondrial genes might also have resulted from 
maternal philopatry and male dispersal from natal groups, 
common in living Trachypithecus and possibly the 
ancestral state of the genus, which could have produced 
these different patterns of genetic inheritance (Karanth et 
al, 2008). 

Previous study including T. p. crepuscula from 
Vietnam and T. p. phayrei from India have found 
Trachypithecus phayrei to be paraphyletic (Karanth et al, 
2008). As a result, it remained unclear whether the 
subspecies T. p. crepuscula is itself monophyletic and 
whether T. p. shanicus is close a relative to either of the 
other two subspecies. For the first time, we found T. p. 
crepuscula from Vietnam, southwestern China, and 
Myanmar are monophyletic and share a young MCRA at 
about 0.41 Ma and T. p. shanicus diverged from T. p. 
phayrei at about 0.37 Ma. However, T. p. shanicus is 
represented by a different and probably more ancient 
PRM1 haplotype than the other two subspecies. This 
might be the result of nuclear introgression from T. p. 
shanicus to the ancestor of T. p. crepuscula due to male 
dispersal and female philopatry. Additional taxon 
sampling of other species and multiple unlinked genes 
are needed for reconstructing these robust evolutionary 
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